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13.  ABSTRACT  (Maximum  200  words) 

We  present  interpolating  functions  which  have  three  orders  of  differentiability  at  each  (convex)  data  point  These  functions 
are  defined  as  piecewise  conics  and  are  therefore  guaranteed  to  be  convex  in  the  case  of  (strictly)  convex  data.  The  modifier 
"affine"  refers  to  the  fact  that  we  make  no  use  of  Euclidean  distance  or  angle  in  the  discussion.  We  also  include  a  discussion 
of  numerical  differentiation  using  conics.  The  nodal  derivatives  for  the  conic  ^lines  satisfy  a  locally  quadiivariate  quadratic 
system  solved  by  Newton  iteration-each  iteration  involving  the  solution  of  a  pentadiagonal  linear  system.  Initial  values  for 
Newton  iteration  are  obtained  by  the  aforementioned  conk  numerical  differentiation.  A  discussion  of  numerical  quadrature 
based  on  conic  splines  is  also  included,  as  well  as  a  discussion  of  what  we  refer  to  as  "sketched"  interpolation,  which  makes  use 
of  the  mathematical  machinery  behind  conic  differentiation  and  local  C?  conic  ^lines.  Sketched  interpolation  is  more  generally 
applicable  than  global  conk  splines  are,  as  well  as  being  computationally  simpler,  more  flexible,  and  smoother  in  a  local 
pointwise  sense.  Hiis  apparent  increase  of  smoothness  beyond  is  obtained  through  a  process  of  le-sketching  during  the 
construction  of  the  interpolant.  Sketched  interpolants  reproduce  conics  with  or  without  re-sketching.  This  is  to  say  that  if  the 
discrete  data  comes  from  a  conk,  all  the  points  of  the  fetched  interpolant  will  lie  on  that  conic. 
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INTRODUCTION 


We  are  mainly  concerned  with  the  smooth  interpolation  of  convex  data.  We  can 
interpolate  nonconvex  data,  but  it  is  obvious  that  we  cannot  get  piecewise  conics  to  have  more 
than  one  continuous  derivative  at  an  inflection  point.  Inflection  points  may  be  few  and  far 
between  anyway.  Twice  differentiable  piecewise  cubics  are  commonly  used  for  interpolation  but 
they  carry  no  guarantee  of  shape  preservation.  They  may  be  forced  to  preserve  shape  in  the 
sense  of  convexity  or  monotonicity,  but  this  may  entail  giving  up  an  order  of  differentiability 
resulting  in  an  interpolant  of  lower  smoothness. 

We  make  use  of  abridged  notation  that  was  popular  among  analytic  geometers  of  the 
nineteenth  century  but  which  seems  subsequently  to  have  fallen  into  disuse  among  analysts  and 
authors  of  modem  texts  on  calculus  and  analytic  geometry.  The  modem  interest  in  practical 
computer-aided  geometric  design  seems  to  have  revived  it  somewhat,  however. 


LINE  MAPS 

The  function 

l(x,y)=  ax*by+c 

where  a,  b,  and  c  are  scalars  and  (^,y)  is  a  vector  in  represents  one  half  of  a  general  two- 
dimensional  affine  transformation.  We  refer  to  it  as  a  line  map  for  short.  If  we  write 

Kx,y>0 

or  simply 

1=0 


we  have  what  is  called  abridged  notation  for  the  equation  of  a  line  in  the  plane.  We  have  only 
one  trivial  theorem  to  prove  regarding  line  maps.  Suppose  vector  P  is  a  real  linear  combination 
of  other  vectors  Q. 


P=Sk,Q, 


where 


Qr(Xpy) 


Hence, 


P  Sk^x^,y^-(^Sk^XpSk^y^ 
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Now  let 


lix,y)=  ax+by+c 


where 


m)=c=m 


Now, 


/(F)  =  aSk^x^+bSk^y^^c 
=  IIk^iax.+by)+c 
=  2k,(KQ)-c)+c 
=  Sk^KQ)Ml-Sk) 


Therefore,  we  have 


/(27ife, .(?,)=  Sk^liQ,)*mil-Sk) 


and  we  see  that  a  line  map  is  a  linear  map  if  Z(0)— 0  and  that  a  line  map  at  least  behaves  in  a 
linear  fashion  if  the  coefficients  of  the  linear  combination  sum  to  unity.  In  any  case,  we  see  how 
a  line  map  operates  on  a  linear  combination  of  vectors. 


CONIC  DERIVATIVES  OR  TANGENTS 

We  can  use  line  maps  to  conveniently  represent  conics  in  abridged  notation  also.  Let 
F„  (i-2  <  n  <  i+2)  represent  a  sequence  of  five  points  in  the  plane  which  appears  to  derive  from 
a  single-valued  convex  function  of  a  single  variable.  From  the  line 

y-yjjj-yi 

X-Xj  Xj-Xj 


we  define  the  line  map 

li^x,y)=(y-yj)(Xj-Xi)-(x-Xi)(yj-yj) 

and 

for  arbitrary  J^J  and  K. 
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Now  select  a  four-element  subsequence  of  the  aforementioned  five-element  sequence  of 
points  and  let  Pj  be  an  element  of  the  subsequence.  Also  let  a  and  A  be  distinct  line  maps  such 
that 


aix,y)= xa^+yOy+OQ,  Aix,y)= xA^+yA^+A^ 


and 


a(Pp=AiPp=0 


If  P  is  a  point  on  a  conic  through  Pj,  we  have 


P-Pj  =(x-xj,y-yp=ix-xp 


(  V— V  \ 


1, 


y-yj 


x-x, 


=  ccv 


Now 

«(P)  =  a{Pj*  <r  v)  =  a{Pp + ara(v)  -  aa(0)  =  a(fl{v)  -a(0)) 


and  similarly 


A(P)=  a(A(v)-A(0)) 


Hence 


If  we  now  let 


r(.P)= 


a(P)  _  Q(v)-Oo 
A{P)~Aiv)-A^ 


P 


we  have 

v-(hyj\riP)-riPp=r^ 

where 

'  A(lo))-Ao  A^^AyyJ 
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Solving  the  last  equation  for  y/,  we  get 


It  is  now  only  necessary  to  compute  for  a  conic  in  order  to  compute  the  derivative  of  or  the 
tangent  to  the  conic  at  Pj. 


We  must  now  define  a  connective  graph  relative  to  the  four  points  of  the  subsequence. 
Define  exactly  four  edges  such  that  each  vertex  has  exactly  two  incident  edges.  Define  a  line 
map  corresponding  to  each  edge  such  that  line  maps  a  and  b  are  associated  with  nonadjacent 
edges  and  line  maps  A  and  B  are  also  associated  with  nonadjacent  edges.  Of  course,  a  and  A 
must  be  associated  with  the  vertex  at  which  we  seek  the  derivative.  The  product  of  line  maps  a 
and  is  a  degenerate  hyperbola  passing  through  the  four  points,  as  is  the  product  of  line  maps  A 
and  B.  Forming  a  linear  combination  of  these  line  map  products  gives  us  "abridged  notation"  for 
a  family  or  pencil  of  non-degenerate  conics  through  these  four  points. 

ab+AAB=0 


Requiring  this  equation  to  hold  for  a  fifth  distinct  point  determines  the  parameter  A  and  the 
unique  conic  through  five  points.  The  point  in  the  original  sequence  but  not  in  the  subsequence 
is  called  P^.  So,  we  compute 

^  a(P,MP,) 

A(P,)B(P,) 

Now,  dividing  the  conic  equation  hy  A,  'fit  have 

rb*XB=Q 


or 


r=- 


b 


Therefore 


XBiPp  a(P^)biP,)B(,Pp 
biPp  ^  A(PpB(PpbiPp 

We  select  the  edges  of  our  connective  graph  so  that  the  denominator  of  /)  is  not  zero.  In 
addition,  we  select  them  so  that  this  is  the  case  even  when  one  of  the  middle  three  points  is 
collinear  with  its  neighbors.  All  the  graphs  reflected  in  Table  1  have  an  hourglass  edge  pattern. 
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Table  1.  Line  Map  Indices 


; 

k 

a 

b 

•  .  ' 

i-2 

i 

1-2,/ +2 

i'-l,i+l 

/-2,/4-l 

1-1,1 4-2 

i-\ 

1-2 

i,i-l-2 

i-l,i-»-2 

1,1 -Hi 

i 

i>2 

i-2,i+l 

i-2,i 

i-l,i+l 

i+1 

/”f-2 

1-1, 1 -1-1 

1-2,1 

1-2,1 +1 

i-l,i 

i+2 

i 

i-2, i+2 

i-l,i>l 

/-I,/ +2 

i-2,i-l-l 

Noting  that 

^x^u~yryj> 


we  have 

I  / 

V*./ 

where  the  Fs  and  Fs  come  from  the  table.  Using  the  table,  we  can  therefore  simply  write  down 
the  formulas  for  differentiating  conics. 

/  _y<-2-y.>2->‘i-2(y<-2-yi*i) 

I  I 

U-Xi*lA-lJ*2A-U*U-2 

y,  yi-i-y.>i->-i-i(yi-i-y.>2) 

•*1-1  “*<+l  “^l-l(*i-l  “*i+2) 

_  _  ^l-l.i+U-2^y+2J-2^W+l,/-l 

'■‘"z  Z  Z 

*i-l.i+2A-2  *y +M-2  *w+2.i-l 
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./_y.-i-yr>‘.(yf-2-yi) 


yr 


I  I  I 

*i--2^,i+2  *i-U+U+2  i-2^*+U 


/  _yi-i-y.>i-''i.i(yi-2-y.>i) 

_^i-lj*U*2^i-2,i,i>2^i-l,U*l 

~  I  I  I 

<-2^+M+2  i-l,y+2  i-2,M+l 


, ,/  _yi-2-y,.2-^i^2(y,-i-y.>2) 

y/+2“' 


^i+2 


h~Xi+l^  ^i-U+1^'  i+2 

^*-1,1+2,!  ^i~2^+l,i  ^i-l,i+l.i+2 


These  formulas  for  numerieal  differentiation  based  on  eonics  can  be  used  only  when  the  data  is 
locally  convex.  One  must  therefore  suffer  the  inconvenience  of  testing  the  data  for  this  property 
before  using  the  formulas.  Further  caution  is  also  needed  in  that  although  these  formulas  yield 
results  superior  to  corresponding  polynomial  based  formulas,  polynomials  at  least  carry  with 
them  the  guarantee  of  single  valuedness,  whereas  conics  are  globally  bivalued.  This  being  the 
case,  one  must  be  sure  that  the  conic  formulas  yield  derivative  values  of  the  proper  sign  at  initial 
and  final  data  points  (j=i-2  or  j=i+2).  In  any  case,  one  can  at  least  appreciate  the  value  of 
abridged  notation  at  this  point. 


PIECEWISE  CONIC  INTERPOLANTS 

Consider  four  consecutive  points  on  a  convex  curve-the  first  and  last  being  fixed.  Define 
line  map  a  relative  to  the  first  two  points,  line  map  b  relative  to  the  last  two  points,  line  map  A 
relative  to  the  second  and  third  points,  and  line  map  B  relative  to  the  first  and  last  points.  A 
family  or  pencil  of  conics  through  these  four  points  is  therefore  given  by 

ab+XAB=0 

Now,  let  the  second  point  approach  the  first  and  the  third  point  approach  the  last.  In  this  case, 
a=0  and  ^>=0  are  lines  tangent  to  the  curve  at  the  two  fixed  points.  Line  mapsy4  and  B  coalesce 
into  a  single  line  map,  say  c.  The  family  of  conics  interpolating  the  curve  is  then  given  by 
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ab+Ac^=0 


Denote  the  two  fixed  data  points  by  and  P.+j,  the  point  of  intersection  of  a=0  and  6=0  by  I, 
and  the  midpoint  of  by  M.  Of  course, 

a(P,)=6(P,,,)=a(i)=6(7)=0 


Define  basis  vectors  u  and  v  by 


Pm-Pi 

v=M-P.=-^—^  -P.=  ‘ 

*2  '2 


An  arbitrary  point  P  on  our  representative  conic  can  therefore  be  expressed  by 

fPy 

P=I+mu+nv=I+m 


Note  that 


1-m 


m-n  m+n 


Therefore, 


Letting 


we  have 


fl(P)= 


m+n 

2 


b(.P)=^b(.P) 


K,=aiP,^,)b(P) 

a(P)b(P)=^(m^-n^ 

4 
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Now  let  point  Q  lie  on  line  c  =  0.  But  since  P,-  and  Pj+i  determine  c  =  0,  we  have 

c(0=O.  <?=(l-r)P,+rPj„ 

Operating  on  Q  with  a  and  b,  we  have 

aiQ)=raiP^,l),  biQ)={l-f)b(P) 


Eliminating  parameter  r,  we  get 


0(0  _m-lKQ) 


or 


aiQ)b(P)^h{Q)a{P^^,)-K^-Q 


We  therefore  define  c(P)  as 

c(P) = a{P)b{P)  +  HP)aiPi^,)- K, 


m+n 

2 


m-n 


K.-K, 

2  '  ' 


From 


a(P)&(F)+^c(P)2=0 

4 


we  have 

— (m^-n^)+  — =0 

4  4 


or  equivalently 

This  is  therefore  the  equation  of  our  family  of  conics  in  local  (m/t)  coordinates.  Note  from  the 
definitions  of  u  and  v  that  we  want  0  <  m  <  1  and  - 1  <  n  <  1.  Also,  Aj/iQ  must  be  negative.  If 
AjA;  <  -1  or  =  -1  or  >  -1,  the  conic  is  an  ellipse,  a  parabola,  or  a  hyperbola,  respectively. 
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If  we  now  differentiate  the  last  equation  three  times  with  respect  to  w,  we  have 

nn=m+XjK.fim-l) 

^+nh=\■^X^^ 

3nn+nn=0 


from  which  we  conclude  that 


nf=3XfKp  nf^^=-3X^f 


In  the  following,  let  the  prime  denote  differentiation  with  respect  to  x  and  the  dot  denote 
differentiation  with  respect  to  m  as  before.  Recalling  that 

P=I+mu+nv 

ac=/j+»iMj+nVj 

y=/2+»i«2+nv2 


we  have  the  following 

y 

/_y_«2+nV2 

X  «j+nvj 
^  («j+jiVj)® 

//-  (^^v^-U2v;m^l*nv;i-3v^n'^) 
dm  («j+nVj)* 

Abbreviating 

ttAv=«iV2-«2Vi 
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we  have  the  first,  second,  and  third  nodal  derivatives  at  the  extremes  of  the  subinterval: 


/  “2-^2 


yi  = 


Uj-Vi 


/ 


«1+Vi 


//_ 


// 

yi*i=- 


(«i+Vir 


(“i-vir 


-3yt^,*Av(«i+Vj+ViA^.) 

3’j.i= - 


(“i+vir 


Now  we  evaluate  the  various  expressions  involving  the  components  of  u  and  v  in  terms  of 
the  data  and  first  derivatives.  We  leave  out  numerous  algebraic  steps  for  the  sake  of  brevity. 
First, 

“2=^2  ■^3’/ (“1  -Vi)  =  -V2  +y4i(“i  +Vi) 

2v2-vi(y/+y4i) 

“1  /  / 


Let 


Hence 


Avi  /  / 

rX-> 

AX; 


Ay/  2Ay/ 
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-v^(A^+B) 

^yl 


-Vl  = 


^y'i 


«1+Vi 


-Vi(A,.+5,)  -i4,Ax 

— - — ■ — -  +v,  = - 

Ay/  Ay/ 


Now 


Ax, 


“2=V2+y/(«i-Vi)= - ^(9iAyi  -25p'j) 

2Ay, 


but 


9,Ay/-2V=-^P'ii-V 


Hence 


“2=- 


-Ax,(y4p?4,+Bp'/) 


2Ay/ 


Also, 


-g,Axf(A,+g,)  Axf(Ap>4i+g,y4_A,g,Ax; 

4Ay/  4Ay/  2Ay/ 


“1^2- “2^1  = 


In  addition. 


I=M-U  =P,+V-K=Pj,.j  -v-« 


Hence 


P.Ax,  ^fAx, 

Ayj  Ay, 


But  also. 


Therefore,  we  have 


Ay,  Ay,- 


We  now  define  a  and  b  explicitly  and  compute 

aix,y)  =(y-y,)(/i  -x)-{x-x)ilj^-y) 


bix,y>iy-yM)(ii-Xt^^ -(x-x^^^)(i2  -yi+i) 


Thus 


Ly^il^-x)- ^xfI^-y)= 


^y'i 


b{P)  =  -  Ay//i  +  Ax,.(/2-y,.,i) = 


^y'i 


-AfsfAx*  ^ 

Kr<P,.,)b{P)-  '  •  <0 

Ay, 

We  see  here  that  A,  must  be  positive.  Now  we  are  prepared  to  evaluate  the  second  and  third 
derivatives  at  the  ends  of  the  subinterval  in  terms  of  the  (x,y)  data,  the  nodal  derivatives,  and  A,. 


yi  =- 


-X^^uKv  1,  .3,  3 


//  _WAv_i  3.  3 

^i+l  ,  ^3  2^'' 


("l+Vj) 


^,„_3A^^Av».-v,-v,W .  V-V 


‘  ^  r«.+vA5  2A..  ^  2 


(“l+Vj) 
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DERIVATIVE  CONTINUITY 


In  this  section,  we  use  the  conditions  for  continuity  of  the  second  and  third  derivatives  to 
derive  a  system  of  equations  solvable  for  the  nodal  derivatives  of  a  piecewise  conic.  Enforcing 
continuity  of  the  second  derivative  at  x,  gives  us 


We  notice  here  that  if  we  were  to  assign  all  the  nodal  derivatives  by  conic  differentiation,  some 
other  local  method,  or  even  arbitrarily,  we  would  only  need  to  assign  A  on  one  subinterval  and 
use  the  previous  continuity  condition  recursively  to  obtain  all  the  other  A’s  for  a  piecewise 
conic.  There  is  definitely  something  to  be  said  for  interpolants  of  lower  smoothness,  especially  in 
the  case  of  piecewise  conics.  For  instance,  we  obviously  do  not  want  any  more  than 
smoothness  to  link  a  curve  with  a  line  in  the  case  of  conics,  since  only  degenerate  conics  can  be 
flat.  But,  since  we  have  advertised  the  ultimate  in  smoothness  for  piecewise  conics,  we  continue 
in  this  direction.  Enforcing  continuity  of  the  third  derivative  at  jc,-  gives  us 


*i-l 


Now,  we  use  the  second  derivative  continuity  condition  to  solve  for  A;.i,  substitute  the  resulting 
expression  into  the  third  derivative  continuity  condition,  simplify,  and  finally  solve  the  resulting 
equation  for  A,. 


X 


i-l 


X/ifAxf 

BliAxl, 


A..jAfB^AxfAXi_iAqi.i 

The  algebra  involved  in  the  last  step  is  mildly  tedious,  but  there  are  numerous  fortuitous 
cancellations.  Now,  enforcing  continuity  of  the  second  and  third  derivatives  at  x.+i  gives  us 


*  "i+l 


In  a  similar  manner,  we  now  eliminate  A^+i,  using  the  second  derivative  continuity  condition,  and 
solve  for  A;  in  the  condition  for  third  derivative  continuity,  getting 


13 


Now,  we  simply  equate  the  two  expressions  defining  A,.,  and  clear  of  fractions,  getting 

-■®i-i^A>iA9iAx,.,i(4(.iAy/Ax,.,-fi,Ay/_iAx^)=0 


This  equation,  which  is  the  main  result  of  this  section,  is  a  necessary  condition  for  continuity  of 
the  second  and  third  derivatives  at  two  neighboring  points  in  terms  of  the  data  and  first 
derivatives  at  four  neighboring  points.  It  behaves  quadratically  in  the  two  inner  derivatives  and 
linearly  in  the  two  outer  derivatives.  We  therefore  refer  to  it  as  a  component  of  a  locally 
quadrivariate  quadratic  system.  It  behooves  us  at  this  point  to  make  a  preliminary  examination 
of  the  system.  We  can  write  a  component  equation  only  for  interior  subintervals,  because  the 
quadrivariate  nature  of  each  component  equation  forces  us  to  exclude  the  first  and  last 
subintervals.  Assuming  that  we  have  n  points  of  data,  we  can  therefore  write  only  n-3 
component  equations.  This  leaves  us  free  to  specify  derivatives  at  the  first  and  last  points  and 
one  interior  point  by  either  conic  differentiation  or  other  means.  If  we  use  conic  differentiation 
to  specify  the  three  derivatives  and  the  data  comes  from  a  conic,  the  conic  is  preserved.  This 
gives  us  a  numerical  check  on  all  the  algebra.  The  need  to  set  one  internal  derivative  causes  the 
linear  systems  involved  in  Newton  iteration  to  be  pentadiagonal. 


SYSTEM  SOLUTION 

This  locally  quadrivariate  quadratic  system  is  solved  using  Newton  iteration.  Initial  values 
for  the  unknown  nodal  derivatives  can  be  computed  using  conic  differentiation.  We  define  the 
following  expressions  for  the  sake  of  brevity. 

Pi  A?  i 

Ay/Ax.vi 


Our  system  component  equations  are  therefore 

fr«iYi-PA=0 

The  system  we  wish  to  solve  is  therefore 

f(y)=o 

augmented  by  three  trivial  equations.  A  Newton  iteration  is  accomplished  by  solving  the 
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linear(ized)  system  (augmented) 

my')-  -fiy') 

for  dy’  and  incrementing  y'  by  dy'.  To  write  this  system,  we  must  now  compute  the  partial 
derivatives  of  the  fs.  Abbreviate 


Compute 


3*-,^r0 


^i-iYrO 

diYr-Bi,Ax.^r^yl*A^i 

di-A=^y'A^i-i*BM 

dA=-KAh-i-BAx. 

a  ,5=0 

^i+2  I 


^r^.BjYi^  YidjCcrfifiAr^A^i 
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The  linear  system  in  which  these  expressions  are  used  is  given  by 

l-^o=0 

9(/irfyo+3j/idy/+a/idy2+a3/i<fy3=-/i 


l<fy/=0 

l-dy'.,=0 


Note  that  in  the  equation,  the  diagonal  element  of  the  matrix  is  associated  with  dy^.j'  and 
there  are  two  nonzero  terms  to  the  right-  Also  note  that  in  the  "Z^"  equation,  the  diagonal 
element  of  the  matrix  is  associated  with  and  there  are  two  nonzero  terms  to  the  left.  The 
'  linear  Newton  iteration  system  is  therefore  pentadiagonal. 


NUMERICAL  QUADRATURE 

We  now  consider  the  use  of  conic  splines  for  numerical  integration  using  previously 
obtained  material.  Considering  a  single  subinterval,  let^l  be  the  area  between  the  conic  and  the 
chord.  Let  P  be  a  point  on  the  conic  and  g  be  a  point  on  the  chord  such  that  P-Q  is  parallel  to 
u.  Our  element  of  area  is  therefore  given  by 

dA^dQKP-Q) 


Now  Q  is  given  by 


Q=Pi*piPi.i-P)=Pi^2pv 


and  P-Q  is  given  by 


P-Q=I+mu+nv-Pf-2pv 
=P.+v-u+mu+nv-P^-2pv 
=(ot-1)m+(1  +n-2p)v 


But  since  P-Q  is  parallel  to  u,  we  have 

u  A(P-0=O=(1  +n-2p)M  Av 


and 


16 


Hence, 


dQ=vdn 
P-Q=(m-\)u 
dA=vdn  Nljn-l)u 
=«Av(l-my&t 


and  iSnally, 


i4=«Av  J  i\-m)dn 
-1 


Defining 


and  recalling  that 


we  can  solve  this  quadratic  equation  for  1-m,  getting 


1 

1-5^ 


If  we  were  now  to  plug  this  expression  into  the  previous  integral  defining  >4,  we  would  be  able  to 
obtain  an  analytic  expression  for  the  integral  in  terms  of  simple  algebraic  and  transcendental 
functions.  Unfortunately,  the  resulting  expressions  are  asymptotically  unrobust  in  the  numerical 
sense.  By  this  we  mean  that  as  the  integration  mesh  becomes  dense  in  the  limit,  s  approaches 
unity  and  produces  numerically  indeterminate  forms.  The  fact  that  s  approaches  unity  in  the 
limit  is  fairly  easy  to  prove,  but  the  algebra  is  somewhat  tedious,  so  we  omit  the  proof  here.  It 
might  be  expected  that  conic  splines  would  become  (oriented)  parabolic  splines  in  the  limit 
anyway.  We  therefore  opt  for  a  quite  robust  infinite  series  development  instead.  Consider  the 
function 


where  x  is  close  to  zero.  Differentiating  successively,  we  have 


2 

(2^-3)L(i-^-¥ 

(jk-2)!  2“"2 

and  therefore 

(t>i) 

t!  *l(*-2)l  2“'® 

We  therefore  have  the  Taylor  series 


ax)=y: 


*=0 

-  /T^-1  (2^-3)!  X  * 

2  fe!(it-2)!  22*^2 


=1 


X  A  (2fc+l)!x*^^ 
2  ife!(jt+2)!  2^**2 


from  which  we  get 


(2A:+1)!  X 


*+2 


2  ife!(ifc+2)!  2^*^ 


Therefore, 


’  2  £S  «(t»2)l2“*“ 


2Nt+2 


and  finally, 

2  h  kKk+2)\  2^"2 


This  series  is  obviously  very  rapidly  convergent  and  robust  in  the  limit  as  s  approaches  unity.  We 
must  now  integrate  this  series  term  by  term  with  respect  to  n. 
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Ul-m)dn 

J  ^  ta  t!(ifc+2)!  2^*2 


=  f  f  (1-n^f^^dn 

i  ^  k!(k+2y.2^*^  i 


Now,  we  must  consider  computing  an  integral  of  the  form 

1 

f  a-x^fdc 
0 


We  do  this  by  sucx^essive  integration  by  parts,  using  the  formula 


b 


b  b  X 


j Jix)g(x)dx=fib)  J  g(x)dx  -  j fix)  f  g(t)dtdx 


to  finally  obtain 


_  2^(p\f 


J  (2p+ 


(2p*iy 


ft’om  which  we  get 


_2^"*((jt+2)!)^ 


(2ife+5)! 

4  2 


f  l-n^dn=—=- 
i  3!  3 


We  therefore  have 


/ 1  ,V  (2**l)Kl-sV‘2“"‘(»^2)l)’ 


which  easily  simplifies  to 


(2/k+3)(21fe+5); 
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But, 


«Av= 


2Lyi 


Hence,  we  finally  have 


A= 


Ay/ 


fE 

^  *=0 


(2jt+3)(2jfe+5)^ 


Now,  if  we  simply  add  in  the  trapezoidal  component  of  the  integral,  we  have  the  following 
general  quadrature  formula 


}<x)ix=— (yj+yj,i) 

*1 


Ay, 


fE 

■5  *=o 


(2ifc+3)(2jfe+5)^ 


Recalling  that 

^  2fl,.^(A,.^Ay/AJ^,.l-B,.Ay/ iAj:,) 

^i-Ai  ^I'A'^j  A'^i-jA^j-.j 
2A,,,(R,,iAy/Ax,.,^-A,Ay4iAx,) 

and  that 


a2d2a  4 

*  ~  A  '2 

Ay, 

we  easily  have  that 

2g,.,g,(i?,Ay/^Ax,-A,.^Ay/Ax,.,) 

^i-iAx^.,A?,..iAy/^ 

^A>i(^iAy4iAx,-g,,^Ay/Ax,,^) 

We  have  three  different  options  with  regard  to  using  this  quadrature  formula.  The  first 
and  most  accurate  option  is  to  compute  the  nodal  derivatives  for  a  conic  spline  and  then  use 
the  quadrature  formula  in  conjunction  with  either  of  the  two  previous  formulas.  The  second  and 
less  accurate  option  is  to  simply  obtain  the  nodal  derivatives  by  conic  differentiation  and  use  the 
average  value  provided  by  the  two  previous  formulas.  The  third  and  least  accurate  option  is 
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again  to  use  conic  differentiation,  but  then  to  simply  discard  the  infinite  series  altogether-- 
especially  if  the  mesh  is  dense.  This  last  option,  though  least  accurate,  is  quite  viable,  since  it 
requires  the  least  computation  and  integrates  any  oriented  parabola  exactly.  In  fact,  these 
quadrature  formulas  are  sufficiently  robust  to  accommodate  infinite  derivatives,  provided  that  the 
numerators  and  denominators  of  the  conic  derivatives  are  calculated  separately  (tangent  vectors). 
Since 


A^y'i-Qr 


dx^  AXf  Ax^dx^  Ax^dx^ 


BryU-qr 


Axt  Ax^dx^^^ 


and 


A  /  /  ! 


dx.  dx^dx^^^ 


^y'i 


AfAdf  ^^2 

Ax^dx^  Ax^dx^^^  '  _id^^A^d^^^^ A) 


The  quadrature  formula  in  this  case  is  easily  seen  to  be 


y(x)dx=—(y^-^yi^i) 


1  (d-AA,)(rf,,iAA,) 
3  dM.i 


It  is  instructive  to  compute 


1 


using  this  quadrature  formula  and 


Ax,.=l,  Ay,=l 

dx^=0,  dXf^i=2, 
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SKETCHED  INTERPOLATION 

By  sketching,  we  mean:  starting  with  at  least  five  data  points  (or  fewer,  if  initial  tangents 
are  specified)  and  during  multiple  passes  over  the  accumulated  data,  computing  the  nodal 
tangents  by  conic  differentiation  and  inserting  exactly  one  point  between  each  two  data  points  on 
each  pass  using  the  previous  expressions  for  XjK^.  We  can  thereby  produce  a  set  of  discrete 
points  as  dense  as  we  wish  to  represent  the  curve.  If  the  original  points  lie  on  a  conic,  the 
generated  points  will  also  lie  exactly  on  that  conic.  In  addition,  this  method  of  sketching 
accommodates  itself  not  only  to  single-valued  data  but  to  double-valued  data  and  closed  curves 
as  well,  because  the  expressions  for  may  be  cast  into  another  more  robust  form.  First,  we 
recall  that 

2g,_tg,(g,Ay/_xAx,-A,.tAy/Ax,..,) 

A,.=(AXj,Ayj),  d^=(dXpdy;),  ia,b)K(c,d)=ad-bc 


Then  we  have 


We  also  have 


»  +  l 


Ax 


i+l 


Ay._  A.AA.,t 

Ax,.  Ax,Ax,.,j 


B,  Ay/.i  Ax.-A,..,  Ay/Ax..i 

ViM-i 


dx^_^dx^  Ax,..i<iJC<.i  dx^dx^,^ 


dx~  *dx  dx.  , 


and 


Aj-iAdj 


i+l 


B 


AXj.ifix. 


Ai-iAA^I 

■‘Ax,..iAxJ 

=dtc.  ^dxdx 

I  I 


(A,.iA^/^(A,Ad,,i) 


i+l 


Then,  multiplying  the  appropriate  quantities,  we  finally  have 
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'  '  (A,.iAd,.i)(A,.iAAMA4*i)' 


2(A,.iAd,)(A,A4,i) 

(A,A4,i)^.iAd,) 

This  is  the  most  robust  form  for  computing  a  value  of  AjiSQ  which  ensures  continuity  of  the  second 
and  third  derivatives  On  the  left.  This  formula  alone  would  be  used  for  the  last  submterval. 

Now  recall  also  that 

.  2.4,^,.,(A,.Ay4tAx,.-g,,iAy/Ax,,t) 

Z.  \  I  *  /  2 


We  then  have  that 

^  AjAdj 

- ((A,A<f^(4,Adi,2)-(AMA^,.2)^A<^j*i)) 

dx^dx^^^dx^^2 


and 


A, A, 


i+1 


AfAdf 

Ax^dXf  AXj^iO^ 


/  2 


^i+1  M,.,  A,AA  i+1  (-iM,,' 

■^^i+r 


Ax,,idXf,2 


AXfAx^ 


■i*l 


dXfdx^ 


'1+1) 


...  (A,ArfJ(A,.,A4i) 

‘  ‘•\A,.,A<i,.^(A,AA,.,X<i,Ad,./ 


Multiplying  the  appropriate  quantities,  we  then  have 
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i,K,. 


(^j+l^^j+2^(^j^^i+l)(^j^^<+l) 

(Aj  A  A  j,i)(4  A^j)  [  (A^,i  Ad,,2)(<^/4,i)  J 


This  is  the  most  robust  form  for  computing  a  value  of  X^,  which  ensures  continuity  of  the  second 
and  third  derivatives  on  the  right.  This  formula  alone  would  be  used  for  the  first  subinterval. 

For  interior  subintervals,  we  compute  a  value  for  A^-  for  the  left  and  right  three  point  local 
conic  splines,  insert  a  point  on  each,  and  take  the  midpoint  of  these.  All  we  need  do  at  this 
point  is  to  show  how  to  compute  a  single  point  on  the  desired  curve.  Letting  P  be  this  point, 
recalling  that 

P=I+mu+nv 

n^=m^+X^^On-lf=m^-s\m-iy' 


and  letting  n  =0  gives  us 


P=J+mu 


But  from 


I=M-u 


we  have  that 


P=M-{\-m)u=M-— 

1+s 

Now,  it  is  only  necessary  to  compute  the  components  of  u.  Recalling  that 

“i""  / 

2Ay/ 

Aj;^(X^/i+gp>/) 

2Ay/ 


we  have  almost  immediately  that 
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2(dMi.t) 


RE-SKETCHING 

The  process  of  simple  sketching  produces  a  dense  set  of  data  quite  smooth  to  the  eye,  but 
an  examination  of  the  second  derivative  or  of  the  curvature  (in  a  Euclidean  context)  reveals  a 
certain  erratic  behavior.  This  microscopically  erratic  behavior  can  be  eliminated,  however,  by  the 
process  of  repeated  re-sketching  during  each  pass  over  the  accumulated  data.  The  process  of  re¬ 
sketching  is  as  follows:  On  a  given  pass  over  the  data  and  after  a  point  has  been  inserted 
between  each  of  the  previous  data  points,  we  recalculate  the  nodal  tangents  at  each  of  the 
previous  data  points  using  the  previous  points  and  the  new  points,  throw  away  the  new  points, 
and  recreate  them  using  the  updated  nodal  tangents  at  the  previous  data  points.  We  can 
obviously  repeat  this  process  of  re-sketching  as  many  times  as  we  wish  per  pass.  Even  only  a  few 
re-sketchings  per  pass  is  sufficient  to  markedly  "iron  out"  the  previous  erratic  behavior  of  the 
second  derivative  or  curvature.  What  results  is  a  pointwise  interpolant  that  api^ars  to  have  even 
greater  smoothness  than  C^.  This  might  be  expected,  however,  since  global  conic  splines  are 
infinitely  smooth  everywhere  except  at  the  original  data  points.  Another  convenient  aspect  of 
sketched  interpolation  is  the  fact  that  we  need  not  specify  any  nodal  tangents  once  and  for  all 
because  of  degree  of  freedom  restrictions,  as  we  must  with  the  global  conic  splines.  On  the 
other  hand,  if  we  wish,  we  can  specify  any  or  all  original  nodal  tangents  once  and  for  all.  The 
inherent  flexibility  of  sketched  interpolation  is  evident.  The  only  thing  that  sketched 
interpolation  does  not  offer  us  is  an  analytic  function  describing  the  interpolant  between  each  of 
the  original  data  points.  But  since  we  often  obtain  such  an  interpolant  only  in  order  to  generate 
a  dense  set  of  points  on  the  interpolant  anyway,  sketching  just  eliminates  the  middleman  and 
gives  us  greater  flexibility  and  smoothness  to  boot.  Since  we  compute  curvature  for  the  sketched 
curves  in  the  Appendix,  we  develop  here  the  formula  for  curvature  at  the  midpoint  of  a  three- 
point  conic  spline  for  reference.  The  general  formula  for  curvature  is  given  by 


Now,  recall  that  the  second  derivative  of  a  conic  spline  at  node  i  is  given  by 
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,,  AfAx^2BaA,_,^y|^x,,,-B,AyU^x,) 

^'2  2 

^A-iiBjAy'-t  Ay/AA:^.,) 

KlBMAx^.^Aq^.^ 

Breaking  this  expression  down,  we  have 

B^AyUAx.-A._iAy'  Ax._^ 

-  Ax 

Ax,dx,,i  *dx,_idx,  Ax..jdx,..i 


- 1  j  + 1 


and 


AjArfj  A,_jAd,. 

_ Ax,.<fa,.  Ax,.t^^ _ 

^*i-i^i-i  Ax,_jAX|. 

‘^i-i«faM  (A,A.f^(A,.,Arfj 

dxf  W,.,A<i,-iXA,A4,Ki,-.Ai^ 


Multiplying  the  appropriate  quantities  together,  we  then  have 


//  1 


Vi 


dXi 


- — ((A,.Ad.,i)(d;  ,Adj)-(A,.  jAd.  MM.i)) 

(A,.iA4.i)(A,A^/,,i)(A,.iAA,)^'  ' 


1  (A.A^ij)(A..,A^^ 

(  d.  ,Ad.  d.Ad. 

1-1  I  1  1+1 

1  1 

A,-,Ai, 

1,  A,-i  A<i._ j  A,. Adj^j  ^ 

CO 

! 

Therefore  the  curvature  is 
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dXi 


1  + 


dx^ 


(dxf  +dyf) 


FINAL  NOTE  ON  TANGENTS 

We  now  present  a  more  straightforward  and  elegant  derivation  of  the  general  formula  for 
a  tangent  to  a  conic.  This  derivation  involves  the  wedge  product,  which  we  have  previously 
defined  as  a  real-valued  binary  operation  on  two-dimensional  vectors: 

mAv  =(Hi,M2)^(Vi,V2)  =  «iV2-a2Vi 

The  wedge  product  is,  of  course,  just  the  oriented  area  of  a  parallelogram  embraced  by  the  two 
vectors.  Note  that  area  and  wedge  product  are  affine  concepts  rather  than  Euclidean  ones. 

Using  the  previous  definition,  it  is  trivial  to  prove  that 

(h  +v)Aw  =«  Aw  +vAw 


and 


d(«  Av)=  duhv  +u^dv 


The  equation  of  a  line  through  two  points  is  given  by 

X-Xf  Xj-Xf 


or 


(y-yi)(Xj-Xf)  -(x-x,)(yj-y;) = 0 
=(Xj-Xpyj-y;)A(x-x,,y-yi) 

=(Pj-P,)A(P-P,) 

=A^A(P-P) 


A  line  map  may  therefore  be  defined  in  terms  of  wedge  product: 

l^iP)=A^A(,P-P) 


fi'om  which  we  immediately  have 
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dl^iP)=A/sdP 


Recall  that  the  equation  of  a  conic  in  abridged  notation  is  • 

a(,P)biPy^MP)BiP)=0 

where  a,  b,A,  and  B  are  line  maps  and  F  is  a  point  on  the  conic.  In  the  case  of  line  map  a,  for 
instance,  we  have 

aiPMPj-Pj)MP-Pt) 

=  a/jHP-Pj)  ‘ 

da(P)=  Kt^P 


Now,  differentiating  the  equation  of  a  conic,  we  have 


a(P)dKPy  KPyiaiP) 
+/l(A(P)dB(P)+  B(P)<i4(P))=0 

a{F)A^MP*KP)A^MP 
*X{A(P)AgSdP^B{P)A^bdP)  =  0 


and  finally 


T(P)KdP=0 


where 

71[P)=a(P)4,+W4,+>l(A(P)4^+B(P)4^) 

Hence,  dP  and  T(P)  are  parallel  and  T{P)  is  a  tangent  vector  to  the  conic.  As  before,  if  P*  is  a 
point  on  the  conic  for  which  none  of  the  line  maps  are  zero,  we  have 

A- 

^(P*)P(P*) 


and 


nP)=a(P)A^+lKP)A^ 
aiP,)b(P,) 


A{P,)BiP,) 


iAiP)A.^B(P)A^) 


Another  tangent  vector  proportional  to  this  one  is  therefore 
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a(P)A,^b{P)A^  A(F)Aj,^BiP)A^ 


In  the  earlier  section  on  conic  differentiation,  we  sought  the  derivative  or  tangent  at  Pj  and 
contrived  to  make 


a(Pp=>4(P;=0 

The  tangent  vector  in  this  case  simplified  essentially  to 

as  it  does  here  also.  However,  if  one  or  two  pairs  of  the  five  original  points  on  the  conic 
coalesce,  however,  which  is  to  say  that  one  or  two  tangents  are  specified  for  a  conic  through  four 
or  three  points  respectively,  then  we  cannot  use  the  simpler  tangent  formula  and  must  resort  to 
the  more  general  one  where  one  or  two  of  the  A’s  would  be  the  specified  tangents  at  the 
coalesced  points.  Note  also  that  before  two  points  are  allowed  to  coalesce,  they  must  be  linked 
by  a  line  map.  Another  important  point  that  should  be  noted  here  is  that  in  the  case  of  strictly 
convex  data,  we  can  use  the  general  formula  for  t(P)  to  calculate  the  tangent  at  any  data  point 
after  having  made  any  reasonable  selection  of  line  maps  and  associated  vectors  once  and  for  all! 
This  eliminates  all  the  special  cases  mentioned  in  the  previous  section  on  conic  differentiation, 
making  the  algebra  more  compact  and  the  programming  more  straightforward. 


CONIC  SECOND  DERIVATIVES 

In  addition  to  obtaining  tangents  to  local  interpolating  conics,  one  may  occasionally  wish 
also  to  obtain  second  derivatives  of  these  interpolants.  Recall  that  the  general  formula  for  a 
tangent  to  an  interpolating  conic  is 


t(P)= 


a(P)A,*b(P)A^  A{P)A^^mp)A^ 

h  U 


where 

4=u(P*)b(P*)  4=>i(P*)B(P^ 


Define  the  component  functionals  for  an  arbitrary  vector  v  as 

x(v) = ^((Vj.v^))  =  Vj  y(v) = y((Vi,Vj)) = 


Since 


t(F)A 

ax 
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we  have,  trivially,  that 


/_y(ri 

~x(t) 


Differentiating 


tKUy)=0 


we  have 


-x(ffdy'=dtht 

ax 


But 


and  since 


*- 

dx 


dx  dx  ‘ 


dA  .  dB  . 

— 4„+— d. 


dx 


dx 


t 

Mt) 


we  have 

Jc(Ox = - j - - } - 

dx  It  h 

Taking  the  wedge  product  of  the  last  equation  with  t,  we  have 

liA/i)iA,Nt)  2(d»(d/r) 
x(t)—At= - ; - - r - 


froni  which  we  finally  obtain  the  general  formula  for  the  second  derivative  of  a  local 
interpolating  conic: 


2 

r(d»(d^Ar)  (d»(d/r)' 

x(jt? 

[  J 
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From  this,  we  immediately  have  a  compact  formula  for  the  magnitude  of  the  Euclidean  curvature 
of  a  local  interpolating  conic: 

*= - ^ - 


DEALING  WITH  INFLECTION  POINTS 

In  order  to  make  this  development  reasonably  complete,  it  is  necessary  to  deal  with  data 
that  may  have  an  occasional  inflection  point.  We  say  that  section  i  of  the  data  has  ain  inflection 
point  if 

WjW2<0 


where 


“"2 

What  we  wish  to  do  in  this  section  is  to  augment  the  data  with  additional  inflection  points  and 
establish  constant  tangents  at  these  inflection  points.  Of  course,  the  resulting  interpolant  is  only 
C*  at  these  points.  The  problem  then  is  twofold.  We  must  establish  a  location  and  a  tangent  for 
each  inflection  point.  For  the  method  we  describe,  we  must  assume  that  the  data  is  convex  for  a 
sufficient  number  of  points  both  preceding  and  following  the  section  in  which  we  wish  to  insert 
the  inflection  point.  Suppose  we  wish  to  insert  an  inflection  point  somewhere  between  P,  and 
Define  the  following  sets. 

|i+l  <x j^i+4) 


Also  let 


L(P)= conic  through 
R(P)=comc  through 


In  addition,  define  line  maps  with  respect  to  and  line  maps  a|j,f>jj,/4s,5jjwith 

respect  to  ir^.  Now,  insert  a  point  U  (as  in  sketching)  on  the  leftmost  section  of  R{P)  and  a 
point  V  on  the  rightmost  section  of  lIP;+i).  If  we  momentarily  consider  the  conic  pairs 
[L(U)Ji{U)]  and  [L(F),i?(F)],  it  becomes  obvious  that  a  unique  point  W  exists  on  the  line 
segment  l/Fsuch  that  L{}V)  and  R{}V)  have  the  same  tangent  (direction)  at  W.  The  point  W  is 
given  by 


31 


W={\-r)U+rV  O^r^l 


and  we  seek  the  value  of  r  for  which 

Ar)=t^(W)Kt^W)=Q 


where 


bjjiW)  ajj^W)  A^{W) 


tJW)= 


bj^iW)  Ai,(W) 


Note  that  we  have  made  use  of  the  general  forniula  for  the  conic  tangent  here  with  W>  The 
zero  off  can  be  found  by  any  number  of  methods.  When  the  proper  value  of  r  has  been 
established,  fFand  t(W)  will  already  have  been  computed. 


REVIEW  AND  CONCLUSIONS 

We  have  shown  how  to  obtain  convexity  preserving  piecewise  conic  interpolating 
functions  having  three  orders  of  differentiability  at  each  data  point.  But,  although  the  title  of 
this  rejxjrt  reflects  this  desired  goal,  the  global  conic  spline  appears  (to  the  author,  at  least)  to 
be  of  somewhat  lesser  value  than  other  elements  of  this  report.  More  important  or  more 
practical  aspects  of  this  work  seem  to  be  the  simple  and  general  formulas  for  conic  numerical 
differentiation,  the  simple  formula  for  oriented  parabolic  numerical  integration  which  can  make 
use  of  conic  numerical  differentiation,  and  sketched  interpolation,  which  makes  use  of  both  conic 
numerical  differentiation  and  the  main  algebraic  machinery  behind  conic  splines.  The  fact 
that  sketched  interpolants  have  apparent  smoothness  in  excess  of  is  due  in  part  to  the  fact 
that  none  of  the  nodal  tangents  (except  the  ones  one  might  wish  to  allow  to  remain  fixed)  are 
finally  established  until  the  last  re-sketch  of  the  last  pass.  Since  sketched  inter^lants  are  not 
analytic  functions,  but  merely  discrete  point  sets,  we  can  say  nothing  more  precise  about  their 
added  smoothness  other  than  to  refer  to  them  as  apparently 

The  author  apologizes  for  the  lack  of  diagrammatic  aids  in  this  report,  but  he  feels  that 
verbal  and  algebraic  precision  more  than  make  up  for  this  lack.  A  picture  is  not  always  worth  a 
thousand  words,  or  even  a  hundred,  if  the  words  are  sufficiently  well-chosen  ones.  For  those 
who  disagree,  however,  we  include  in  the  Appendix  some  plots  of  sketched  interpolants  with 
various  levels  of  re-sketching  and  accompanying  polar  plots  of  curvature  with  respect  to  the 
centroid  of  the  region  enclosed  by  the  interpolant. 
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APPENDIX 
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Original  Data  and  Sketched  Interpolant 
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Sketched  Interpol 
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Sketched  Interpola 
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Sketched  Interpol 
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